Critical Scaling of Shear Viscosity at the Jamming Transition 
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We carry out numerical simulations to study transport behavior about the jamming transition of 
a model granular material in two dimensions at zero temperature. Shear viscosity n is computed as 
a function of particle volume density p and applied shear stress a, for diffusively moving particles 
with a soft core interaction. We find an excellent scaling collapse of our data as a function of the 
scaling variable cr/|p c — p\ A , where p c is the critical density at a — ("point J"), and A is the 
crossover scaling critical exponent. We define a correlation length £ from velocity correlations in 
the driven steady state, and show that it diverges at point J. Our results support the assertion that 
jamming is a true second order critical phenomenon. 
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In granular materials, or other spatially disordered 
systems such as colloidal glasses, gels, and foams, in 
which thermal fluctuations are believed to be negligible, 
a jamming transition has been proposed: upon increas- 
ing the volume density (or "packing fraction" ) of parti- 
cles p above a critical p c , the sudden appearance of a 
finite shear stiffness signals a transition between flowing 
liquid and rigid (but disordered) solid states It has 
further been proposed by Liu and Nagel and co-workers 
0, E| that this jamming transition is a special second or- 
der critical point ( "point J" ) in a wider phase diagram 
whose axes are volume density p, temperature T, and ap- 
plied shear stress a (the latter parameter taking one out 
of equilibrium to non-equilibrium driven steady states). 
A surface in this three dimensional parameter space then 
separates jammed from flowing states, and the intersec- 
tion of this surface with the equilibrium p — T plane at 
a = is related to the structural glass transition. 

. Ilfjl. t heoretical 
lEllZlGj] works 
have investigated the jamming transition, mostly by con- 
sidering behavior as the transition is approached from the 
jammed side. In this work we consider the flowing state, 
computing the shear viscosity rj under applied uniform 
shear stress. Previous works have simulated the flow- 
ing response to applied shear in glassy systems at finite 
temperature 3, 20, 21 1, and in foams [j] and granular 
systems [l(| at T = 0, p > p c . Here we consider the 
p — a plane at T = 0, showing for the first time that, 
near point J, rj~ 1 (p,a) collapses to a universal scaling 
function of the variable o~/\p c — p\ A for both p < p c and 
p > p c . We further define a correlation length £ from 
steady state velocity correlations, and show that it di- 
verges at point J. Our results support that jamming is a 
true second order critical phenomenon. 

Following O'Hern et al. Q , we simulate frictionless soft 
disks in two dimensions (2D) using a bidisperse mixture 
with equal numbers of disks of two different radii. The 



radii ratio is 1.4 and the interaction between the particles 
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where is the distance between the centers of two par- 
ticles i and j, and dij is the sum of their radii. Particles 
are non-interacting when they do not touch, and inter- 
act with a harmonic repulsion when they overlap. We 
measure length in units such that the smaller diameter 
is unity, and energy in units such that e = 1. A system 
of N disks in an area L x x L y thus has a volume density 



A%(0.5 2 + 0.7 2 )/(2L x L y ) 
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To model an applied uniform shear stress, a, we first 
use Lees-Edwards boundary conditions to introduce 
a uniform shear strain, 7. Defining particle z's position as 
Vi = {xi + jyi, yi), we apply periodic boundary conditions 
on the coordinates Xi and in an L x x L y system. In 
this way, each particle upon mapping back to itself under 
the periodic boundary condition in the y direction, has 
displaced a distance Ax = ^Ly in the x direction, re- 
sulting in a shear strain Ax/L y = 7. When particles do 
not touch, and hence all mutual forces vanish, xi and yi 
are constant and a time dependent strain "/(t) produces a 
uniform shear flow, dvi/dt = yi(dy /dt)x. When particles 
touch, we assume a diffusive response to the inter-particle 
forces, as would be appropriate if the particles were im- 
mersed in a highly viscous liquid or resting upon a rough 
surface with high friction. This results in the following 
equation of motion, which was first proposed as a model 
for sheared foams 
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The strain 7 is then treated as a dynamical variable, 
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obeying the equation of motion, 
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dt 7 
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where the applied stress a acts like an external force on 
7 and the interaction terms V(r%j) depend on 7 via the 
particle separations, iy, = ([xi-Xj] L:i .+~/[yi-yj] Ly ,[yi- 
yj]z y ), where by [. . .]l we mean that the difference is 
to be taken, invoking periodic boundary conditions, so 
that the result lies in the interval (—L ll /2,L li /2]. The 
constants D and D 1 are set by the dissipation of the 
medium in which the particles are embedded; we take 
units of time such that D = _D 7 = 1. 

In a flowing state at finite a > 0, the sum of the inter- 
action terms is of order O(N) so that the right hand side 
of Eq. (j4]) is O(l). The strain 7(4) increases linearly in 
time on average, leading to a sheared flow of the particles 
with average velocity gradient dv x /dy = (dj/dt), where 
v x {y) is the average velocity in the x direction of the par- 
ticles at height y. We then measure the shear viscosity, 
defined by, 



V 



dv x /dy (dj/dt) 



(5) 



We expect r/" 1 to vanish in a jammed state. 

We integrate the equations of motion, Eqs. 
starting from an initial random configuration, using the 
Heuns method. The time step At is varied according to 
system size to ensure our results are independent of At. 
We consider a fixed number of particles N, in a square 
system L = L x — L y , and vary the volume density p by 
adjusting the length L according to Eq. |(5J). We simulate 
for times itot such that the total relative displacement per 
unit length transverse to the direction of motion is typ- 
ically 7(ftot) ~ 10, with 7(ttot) ranging between 1 and 
200 depending on the particular system parameters. 

In Fig. [TJ we show our results for 77 _1 using a fixed 
small shear stress, a — I0 -5 , representative of the a — > 
limit. Our raw results are shown in Fig. [T^t for several 
different numbers of particles N from 64 to 1024. Com- 
paring the curves for different N as p increases, we see 
that they overlap for some range of p, before each drops 
discontinuously into a jammed state. As N increases, the 
onset value of p for jamming increases to a limiting value 
p c ~ 0.84 (consistent with the value for random close 
packing and r/^ 1 vanishes continuously. For finite 
N, systems jam below p c because there is always a fi- 
nite probability to find a configuration with a force chain 
spanning the width of the system, thus causing it to jam; 
and at T = 0, once a system jams, it remains jammed 
for all further time. As the system evolves dynamically 
with increasing simulation time, it explores an increasing 
region of configuration space, and ultimately finds a con- 
figuration that causes it to jam. The statistical weight 
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FIG. 1: (color online) a) Plot of inverse shear viscosity 77 -1 
vs volume density p for several different numbers of particles 
N, at constant small applied shear stress a — 10 -5 . As TV 
increases, one see jamming at a limiting value of the density 
p c ~ 0.84. b) Log-log replot of the data of (a) as r/^ 1 vs 
p c — p, with p c = 0.8415. The dashed line has slope /3 — 1.65 
indicating the continuous algebraic vanishing of r] -1 at p c with 
a critical exponent /3. 



of such jamming configurations decreases, and hence the 
average time required to jam increases, as one either de- 
creases p, or increases N In the limit /V — > 00, we 
expect jamming will occur in finite time only for p > p c . 
In Fig. [lb we show a log-log plot of 77 -1 vs p c — p, us- 
ing a value p c = 0.8415. We see that the data in the 
unjammed state is well approximated by a straight line 
of slope (3 = 1.65, giving r\~ x ~ \p — p c \^ in agreement 
with the expectation that point J is a second order phase 
transition. 

If point J is indeed a true critical point, one expects 
that its influence will be felt also at finite values of the 
stress cr, with obeying a typical scaling law, 



V- 1 (P,°) = \P-Pcff±( ] ° , A ) 

\|p-Pc| A y 



(6) 



Here z = a/\p—p c \ A is the crossover scaling variable, A is 
the crossover scaling critical exponent, and f-(z), f+(z) 
are the two branches of the crossover scaling function for 
p < p c and p > p c respectively. 

In Fig. [5] we show a log-log plot of inverse shear vis- 
cosity r\~ x vs applied shear stress cr, for several different 
values of volume density p. Our results are for systems 
large enough that we believe finite size effects are negli- 
gible. We use N = 1024 for p < 0.844 and N = 2048 for 
p > 0.844. Again we see that p c ~ 0.8415 separates two 
limits of behavior. For p < p Cl log?? -1 is convex in logc, 
decreasing to a finite value as a — > 0. For p > p c , log?? -1 
is concave in log cr, decreasing towards zero as a — » 0. 
The dashed straight line, separating the two regions of 
behavior, indicates the power law dependence that is ex- 
pected exactly at p = p c (see below). Similar power 
law behavior at p c was recently found in simulations of a 
three dimensional granular material [231 ] - 
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FIG. 2: (color online) Plot of inverse shear viscosity r\~ x vs 
applied shear stress a for several different values of the vol- 
ume density p. The dashed line represents the power law 
dependence expected exactly at p — p c and has a slope 
13/ A = 1.375. Solid lines are guides to the eye. Points la- 
beled a = 0.0012 correspond to densities p = 0.870, 0.872, 
0.874, 0.876, and 0.878. 
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FIG. 3: (color online) Plot of scaled inverse viscosity 7/ _1 /|p — 
p c |' 3 vs scaled shear stress z = a j\p — p c \ A for the data of 
Fig. [2] We find an excellent collapse to the scaling form of 
Eq. © using values p c = 0.8415, /3 = 1.65 and A = 1.2. 
The dashed line represents the large z asymptotic dependence, 
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Data point symbols correspond to those used in 



In Fig. [3J we replot the data of Fig. [2] in the scaled 
variables i]~ 1 /\p — ( o c |' 3 vs cr/\p — p c \ A . Using p c = 0.8415, 
(3 = 1.65 (the same values used in Fig. QJ)) and A = 1.2, 
we find an excellent scaling collapse in agreement with 
the prediction of Eq. ([6]). As the scaling variable z — + 0, 



f—(z) — * constant; this gives the vanishing of rj 



P< 



V 3 at a = 0. As 
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oo, both branches of the scaling 



function approach a common curve, f±{z) 
that precisely at p = p c , r/^ 1 



as a 
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This is shown as the dashed line in both Figs. [3] and [U A 
similar scaling collapse of r\ has been found in simulations 
[2oT | of a sheared Lennard- Jones glass, as a function of 
temperature and applied shear strain rate 7, but only 
above the glass transition, T > T c . By comparing the 
goodness of the scaling collapse as parameters are varied, 
we estimate the accuracy of the critical exponents to be 
roughly, /3 = 1.7 ± 0.2 and A = 1.2 ± 0.2. 

That the crossover scaling exponent A > 0, implies 
that a is a relevant variable in the renormalization group 
sense, and that critical behavior at finite a should be 
in a different universality class from the jamming tran- 
sition at point J (i.e. a = 0). The nature of jamming 
at finite a > will be determined by the behavior of 
the branch of the crossover scaling function f+(z), that 
describes behavior for p > p c . From Fig. [3J we see that 
f+{z) is a decreasing function of z. If f+(z) vanishes only 
when z — > 0, then Eq. © implies that 77 -1 vanishes for 
p > p c only when a = 0, and so there will be no jamming 
at finite a > 0. If, however, f+(z) vanishes at some fi- 
nite zq, then rf 1 will vanish whenever a/ (p — p c ) A — zq; 
there will then be a line of jamming transitions emanat- 
ing from point J in the p — a plane given by the curve 
p*(cr) = p c + (<t/zq) 1 / a . If f+(z) vanishes continuously 



at zq, jamming at finite a will be like a second order 
transition; if f+(z) jumps discontinuously to zero at Zo, 
it will be like a first order transition. Such a first order 
like transition has been reported in simulations 2^, 21 1 
of sheared glasses at finite temperature below the glass 
transition, T < T c . However, recent simulations [101 ] of 
a granular system at T — 0, p > p Cl showed that a sim- 
ilar first order like behavior was a finite size effect that 
vanished in the thermodynamic limit. With these obser- 
vations, we leave the question of criticality at finite a to 
future work 

The critical scaling found in Fig. [3] strongly suggests 
that point J is indeed a true second order phase transi- 
tion, and thus implies that there ought to be a diverg- 
ing correlation length £ at this point. Measurements 
of dynamic (time dependent) susceptibilities have been 
used to argue for a divergent len gth scale in both the 
thermally driven glass transition [25j, and the density 
driven jamming transition Here we consider the 

equal time transverse velocity correlation function in the 
shear driven steady state, 



g(x) = (vy(x l ,y l )v y (x i + x,yi)} 



(7) 



where v y (xi,yi) is the instantaneous velocity in the y 
direction, transverse to the direction of the average shear 
flow, for a particle at position (xi, j/j). The average is over 
particle positions and time. In the inset to Fig. |4]we plot 
g(x)/g(0) vs x for three different values of p at fixed a = 
10~ 4 and number of particles N — 1024. We see that g{x) 
decreases to negative values at a well defined minimum, 
before decaying to zero as x increases. We define £ to be 
the position of this minimum. That g(£) < 0, indicates 
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FIG. 4: (color online) Inset: Normalized transverse velocity 
correlation function g(x)/g(0) vs longitudinal position x for 
N = 1024 particles, applied shear stress a — 10 -4 , and vol- 
ume densities p — 0.830, 0.834 and 0.838. The position of 
the minimum determines the correlation length £. Main fig- 



ure: Plot of scaled inverse correlation length £ 
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scaled shear stress z = <r/\p — p c | for the data of Fig. [21 We 
find a good scaling collapse using values p c = 0.8415, A = 1.2 
(the same as in Fig. [3j and v — 0.6. Data point symbols 
correspond to those used in Fig. [3] 



that regions separated by a distance £ are anti- correlated. 
We can thus interpret the sheared flow in the unjammcd 
state as due to the rotation of correlated regions of length 
£. Similar behavior, leading to a similar definition of 
£, has previously been found (2(| in correlations of the 
nonaffine displacements of particles in a Lennard-Jones 
glass, in response to small elastic distortions. 

As with viscosity, we expect the correlation length 
£(p, a) to obey a scaling equation similar to Eq. ©. We 
consider here the inverse correlation length which 
like ?7 _1 should vanish at the jamming transition, obey- 
ing the scaling equation, 



6 1 (p,°~) = \P- Pc\"hi 



(8) 



The correlation length critical exponent is is, but the 
crossover exponent A remains the same as for the vis- 
cosity. 

In Fig. H] we plot the scaled inverse correlation length, 
£~V|p — p c \ v vs the scaled stress, cr/\p — p c \ A - Using p c = 
0.8415 and A = 1.2, as was found for the scaling of 
we now find a good scaling collapse for £ _1 by taking the 
value v = 0.6. By comparing the goodness of the collapse 
as v is varied, we estimate v — 0.6±0.1. From the scaling 
equation Eq. ([8]) we expect both branches of the scaling 
function to approach the power law h±(z) ~ z v ' A as 
z — > oo, so that c; -1 ~ tW A as a — > at p = p c [23 |. 
This is shown as the dashed line in Fig. [4j Our result 
is consistent with the conclusion u v is between 0.6 and 
0.7" of Drocco et al. 0] for the flowing phase, p < p c . 



It also agrees with v = 0.71 ± 0.08 found by O'Hern et 
al. Q from a finite size scaling argument. Wyart et al. 
14j have proposed a diverging length scale with exponent 
v = 0.5 by considering the vibrational spectrum of soft 
modes approaching point Jfrom the jammed side, p > p c . 
While our results cannot rule out v = 0.5, our scaling 
collapse in Fig. [4] does seem somewhat better when using 
the larger value 0.6. 
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